## Import data from .RData
load("data/cont_snotel.RData")

# SNOTEL Analysis, calculate the number and amount of positive SWE between 
cont_snotel_spring <- cont_snow_data %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
  mutate(newsnow = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
  filter(state != "SD", year >= 1979, ifelse(leap_year(year) == FALSE, yday >= 60 & yday <= 243, yday >= 61 & yday <= 244))

# Determine average number of days with increased SWE per spring
cont_snotel_site <- cont_snotel_spring %>%
  group_by(site_id, site_name, state, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), years = length(unique(year)), avg_spring_days = newsnow_days/years, newsnow_depth = sum(newsnow, na.rm = TRUE), avg_spring_snow = newsnow_depth/years, norm_spring_snow = avg_spring_snow/911.8571)
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
mapview(cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "avg_spring_days", crs = 4326, grid = FALSE)
mapview(cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "avg_spring_snow", crs = 4326, grid = FALSE)
cont_snotel_site_yr <- cont_snotel_spring %>%
  group_by(site_id, site_name, state, year, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), newsnow_depth = sum(newsnow, na.rm = TRUE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
spring_snow_days <- ggplot(cont_snotel_site_yr, aes(x = state, y= newsnow_days, color = state)) +
  geom_boxplot()+
  ggtitle("Days of increased SWE at all sites (1978-2022)")

ggplotly(spring_snow_days)
spring_snow_depth <- ggplot(cont_snotel_site_yr, aes(x = state, y= newsnow_depth, color = state)) +
  geom_boxplot()+
  ggtitle("")

ggplotly(spring_snow_depth)
# spring_snow_days_yr <- ggplot(cont_snotel_site_yr, aes(x = year, y = newsnow_days, group = state, color = state)) + 
#   geom_line()
# 
# ggplotly(spring_snow_days_yr)